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WAVELET  TRANSFORM  FOR  TIME-FREQUENCY 
ANALYSIS  OF  THE  VIBRATIONAL  SIGNATURE 
AND  ITS  APPLICATION 

By 

Jae-Jin  Jeon  and  Young  S.  Shin 

ABSTRACT 

Wavelet  transform  is  applied  to  the  analysis  of  vibration 
signatures  in  order  to  verify  the  ability  of  the  detection  of 
abnormal  condition.  It  can  well  describe  the  dynamics  of  the 
signal's  spectral  composition  of  a  non-stationary  and 
stationary  signal  to  be  measured  and  presented  in  the  form  of 
3-D  time-frequency  map.  Although  wavelet  has  been  developed 
over  about  ten  years  in  the  mathematics  and  physics,  its 
engineering  applications  is  a  first  stage.  The  objective  of  this 
report  outlines  the  definition  of  the  wavelet  transform  and  is  to 
discuss  the  properties  of  the  wavelet  transform  as  new  tool  for 
the  vibration  analysis,  and  then  demonstrates  how  it  may  be 
applied  to  the  machinery  condition  monitoring. 


I.  INTRODUCTION 

Wavelets  are  a  very  popular  topic  of  the  signal  processing  and  applied 
mathematics.  In  the  last  ten  years,  an  interest  in  them  has  grown  at  an 
fast  rate  in  signal  and  numerical  analysis  [Beylkin,  Coifinan  and  Rokhlin, 
1991,  Heil,  1990  and  Resnikoff  and  Burrus,  1990].  Wavelet  analysis  appears 
to  new  subject  for  the  time-frequency  analysis  of  the  vibrational  signatures. 
Tra'^itional  spectral  analysis  provides  spectral  values  which  are 
independent  of  time.  It  is  assumed  to  be  ergodic  and  stationary  signal  with 
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time.  However  the  signal  associated  with  most  vibrational  phenomena  are 
in  general  time  varying,  which  means  that  their  characteristics  change 
with  time  and  they  have  various  features  in  different  time  frames.  For 
example,  the  vibration  during  the  start-up  of  an  engine  or  pump  is  non¬ 
stationary,  the  sound  pressure  generated  from  speaker  is  nonstationary, 
and  so  on.  In  case  of  the  signal  containing  some  transient  or  nonstationary 
conditions,  the  traditional  approach  in  signal  analysis  fails  to  describe  the 
dynamics  of  the  signal's  frequency  components  changes.  Changes  in  the 
condition  of  a  component  such  as  a  gear  can  be  expected  to  cause  some 
change  in  the  vibration  generated  by  mechanical  system.  Very  little 
damage  detection  can  be  performed  using  the  vibration  signal  directly  from 
the  system  because  the  small  changes  generated  by  early  damage  may  be 
masked  by  the  normal  vibration  of  the  system. 

Two  methods  of  signal  analysis  for  nonstationary  application  are 
commonly  used  in  time-frequency  domain.  The  windowed  Fourier 
transform,  otherwise  so  called  short  time  Fourier  transform(STPT),  has  a 
short  time  window  of  a  fixed  size  centered  at  time  t  as  figure  1(a).  The 
windowed  Fourier  transform  is  given  as  follows 

F(oj,t)  =  Jg(T  -  t)  e'‘®^s(T)  dr  (1) 

where  F(o),t)  is  the  windowed  Fourier  transform,  g(t)  window  fxmction  and 
s(t)  time  signal.  The  range  of  integration  is  from  -oo  to  oo.  If  the  length  of 
the  window  is  time  duration  T,  the  resolution  of  time  frequency  domain 
depends  on  T.  Its  frequency  bandwidth  or  frequency  resolution  is 
approximately  l/T.  Therefore  this  method  have  the  limitation  of  resolution 
in  both  time  and  frequency  domain  simultaneously. 

The  second  method,  often  called  Wigner  distribution  method,  is  based 
on  the  instantaneous  power  spectra  defined  as  following  equation  (2) 

w{(o,t)  =  —  j5(/- t/2)  r/2)  dr  (2) 

2n  •' 
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where  w(cu,t)  is  Wigner  distribution  function.  Equation  (2)  is  called  the 
Wigner  distribution  of  5(/)in  figure  Kb).  We  discussed  well  at  references 
[Jeon  and  Shin  (1993)  and  Shin,  Jeon  and  Spooner  (1993))  about  the 
characteristics  of  Wigner-Ville  distribution.  For  a  nostationary  signal 
analysis,  spectrogram  is  commonly  used,  which  is  based  on  the 
assumption  that  it  is  a  collection  of  a  short  duration  stationary  signal. 
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Figure  1.  Calculation  of  the  windowed  Fourier  transform  and 
instantaneous  spectral  density(Wigner  distribution)  of  a  vibration  signal. 


A  major  drawback  of  this  approach  is  that  the  frequency  resolution  is 
directly  affected  by  the  duration  of  short  stationary  time,  which 
subsequently  determines  the  time  resolution.  The  frequency  and  time 
resolution  of  the  Wigner  distribution  are  not  determined  by  the  short 
diiration  but  rather  determined  by  the  selection  of  desired  resolution  of  the 
signal  itself,  but  may  not  be  appropriate  for  signals  containing  patterns  at 
both  large  and  very  small  scales. 

The  windowed  Fourier  transform  is  well  portray  the  characteristics  of 
signals  in  which  all  of  the  patterns  appear  at  approximately  the  same 
scale,  but  may  not  be  appropriate  for  signals  containing  patterns  at  both 
large  and  small  scales  because  of  the  time  fixed  window  size.  The 
multifrequency  channel  decomposition,  which  are  an  intermediate  between 
a  time  and  a  Fourier  representation,  have  found  many  applications  in 
signal  and  image  processing  [Mallat,  1989].  Much  recent  research  has 
been  focused  on  this  domain  with  the  modeling  of  a  new  decomposition 
called  the  wavelet  transform,  which  is  presented  signal  by  summation  of 
family  of  functions  which  are  the  dilations  and  translations  of  a  imique 
function  called  a  wavelet.  Instead  of  portraying  a  signal  into  harmonic 
functions  in  Fourier  transform),  the  signal  is  presented  into  a  series  of 
orthogonal  basis  functions  of  finite  length.  Each  wavelet  is  located  at  a 
different  position  on  the  time  axis.  At  the  finest  scale,  wavelets  may  be  very 
short  indeed;  at  a  coarse  scale,  they  may  be  very  long.  Alternatively  very 
small  disturbances  in  a  record  of  machinery  vibration  can  be  easily 
characterized  from  a  wavelet  S-D  or  2-D  map  in  which  the  mean-square 
value  of  the  time  record  is  shown  over  wavelet  scale  and  position. 

One  important  property  of  a  wavelet  transform  is  its  ability  to 
characterize  easily  the  local  regularity  of  a  function.  Simply  by  a  change  of 
the  scale  parameter(dilation)  in  the  wavelet  transform,  many  scales  of  local 
structure  can  be  described  by  a  distribution  in  the  time-scale  plane.  Wang 
and  McFadden  (1993)  introduced  the  advantage  for  examining  the  vibration 
signal  generated  by  a  gear  and  D.  E.  Newland  (1993)  well  investigated  the 
properties  of  the  wavelet  as  a  new  tool  for  the  analysis  of  vibration  records. 
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This  report  outlines  the  definition  of  the  wavelet  transform,  compare 
with  the  advantages  and  the  disadvantages  of  wavelet  and  pseudo  Wigner- 
Ville  distribution  in  time-frequency  domain  analysis,  and  then 
demonstrates  how  it  may  be  applied  to  analysis  of  the  vibration  signals  for 
the  machinery  condition  monitoring. 
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II.  DEFINITION  OF  WAVELET  TRANSFORM 


To  overcome  the  limitation  of  the  fixed  resolution  of  the  windowed 
Fourier  transform  in  the  time  and  frequency  domains  by  decomposing  the 
signal  s(t)  into  a  family  of  functions  which  are  the  translation  and  dilation 
of  an  unique  function  y/it),  defined  the  continuous  wavelet  transform  as 


W(a,b) 


»od  • 

J  )  s(t)  dt 

a 


(3) 


where  W(a,b)  is  wavelet  transform,  \f/  is  an  analyzing  wavelet,  a 
represents  a  time  dilation,  b  a  time  translation,  and  bar  for  complex 
conjugate.  The  normalization  factor  l/-\/a  is  perhaps  most  effectively 
visualized  as  endowing  |W(a,h)|  with  unit  of  power/Hertz  [Shensa,  1992]. 


We  consider  the  space  L2(R)  of  measurable  function  y/,  defined  on  the 

real  line  R  (R  (-€»,  »>)),  certain  weak  'admissibility'  conditions  are  usually 
required  on  y/(.t)  [Shensa,  1992]; 


i  kl 


6(0  < 


OO 


(4) 


where  ^(<y)is  the  Fourier  transform  of  ^^(t).  They  enstire  that  the 

transformation  is  a  bounded  invertible  operator  in  appropriate  space 
[Daubechies,  1988].  If  yf(Q))  is  differentiable,  then  it  suffices  that  y/  be 

zero  mean,  that  is, 


J  r«)  dt  =  0 


(5) 


for  equation  (4)  to  be  satisfied.  In  particular,  since  the  local  average  values 
of  every  function  in  L2(R)  must  'decay'  to  zero  at  ±«>,  the  sinusoidal(wave) 
functions  (basis  of  Fourier  transform)  do  not  belong  to  l2(R).  Fourier 

series  representation  of  any  functions  is  in  L2(0,2r).  In  fact,  if  it  looks  for 
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wavelets  that  generate  L2(R),  these  wavelets  should  decay  to  zero  at  ±«>  ; 
and  for  all  practical  purposes,  the  decay  shouJd  be  very  fast  [Chui,  1992]. 

In  signal  processing,  the  significance  of  equation  (3)  is  well  understood 
by  comparing  it  to  the  windowed  Fourier  transform  (or  short-time  Fourier 
transform): 


F((o,b)  =  Jg(t  -  h)e‘‘“s(t)dt.  (6) 

oe 

Thus,  to  obtain  F(0),  b),  one  multiplies  the  signal  by  an  appropriate  window 
g(t)  (such  as  Gaussian)  centered  at  time  b  and  then  takes  the  Fourier 
transform.  In  mathematical  forms,  equation  (6)  is  an  expansion  of  the 
signal  in  terms  of  family  of  functions  g(t-b)e'“’,  which  are  generated  from  a 
single  function  g(t)  through  translations  b  in  time  and  translations  0)  in 
frequency.  In  contrast,  the  wavelet  transform  of  equation  (3)  is  an 
expansion  in  function  V'((t  -b)!  a)  generated  by  translation  b  and  dilation 

a  in  time.  Thus,  the  continuous  wavelet  transform  is  similar  to  windowed 
Fourier  transform  with  a  different  window  size  for  each  frequency.  The 
important  facts  of  this  is  that,  while  the  basis  functions  in  equation  (6)  have 
the  same  time  and  frequency  resolution  at  all  points  of  the  transform  plane, 
those  of  wavelet  transform  have  the  time  resolution  which  decrease  with 
increasing  a  and  the  frequency  resolution  which  increase  with  decreasing 
a  width  adapted  to  their  frequency  components:  at  high  frequency  \ff  are 
very  narrow,  while  at  low  frequency  y/  are  much  broader.  As  a  results,  the 
wavelet  transform  is  better  than  the  windowed  Fourier  transform  to 
analyze  on  very  small  dist\irbance,  i.e.,  high  frequency  phenomena.  This 
property  can  be  a  best  advantage  in  signal  processing  since  high  frequency 
characteristics  are  generally  highly  localized  in  time  whereas  slowly 
var)ring  signals  require  good  low  frequency  resolution  [Shensa,  1992]. 
Figure  2  shows  the  typical  shapes  of  windowed  Foxirier  transform  functions 
and  wavelet. 

The  wavelet  transform  means  that  signal  s(t)  is  characterized  by 
decomposition  into  a  set  of  wavelet  family  with  series  of  different  frequency 
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bandwidth.  For  the  decomposition  of  the  signal  s(t)  defined  on  real  line,  it 
is  necessary  to  shift  l//  along  R.  Let  Z  denote  the  set  of  integers; 


Z  =  { . .1,  0,  1. . }. 

The  simplest  way  for  y/  to  cover  all  of  R  is  to  consider  all  the  integral  shift 
of  y/. 


y/(t  -  k),  k  e  Z. 


(7) 


Figure  2.  The  typical  shapes  of  (a)  windowed  Fourier  transform  function 
g(t),  (b)  wavelet  V^(t)  =  (l-t^)e‘*  ^^[Daubechies,  19923. 
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Next,  as  in  Fourier  transform,  it  must  be  also  considered  wavelets  with 
different  frequencies.  It  is  considered  wavelets  with  frequencies  partitioned 
into  consecutive  'octaves'  (or  frequency  bands).  For  computational 
efficiency,  it  will  be  used  integral  power  of  2  for  frequency  partitioning,  that 
is,  the  small  wavelets  is  considered  as  follows 

-  k),  j,k  €  Z.  (8) 

From  the  equation  (3),  \f/(2h  -  k)  is  obtained  from  the  wavelet  function 
y/ii)  by  dilation  of  1/2^  and  translation  of  VJlK  An  important  particular 
case  of  the  discrete  wavelet  transform  which  was  found  is  that  some 
wavelets  y/{t)  exist  such  that  V?  “  2~^k))(j  orthonormal 

basis  of  L2(R)  [Mallat,  1989]. 

As  originally  proposed  by  Morlet  et  al.  (1982),  \f/  was  a  modulated 
Gaussian 


^(t)  e‘“o'  (9) 

and  this  function  is  selected  to  analyr*  5  the  vibration  of  gear  box  for  signal 
processing  applications[Wang  and  McFadden,  1993].  The  example  of  a 
shape  for  the  modulated  Gaussian  wavelet  is  shown  in  figiire  3.  The 
Fourier  transform  of  the  first  wavelet  family  Ij/it/a),  i.e.,  no  translation,  is 


A 

y/iaco)  -  at 


~((0~  (Oo/afa^/2 


(10) 


which  has  analysis  frequency  (OqIq.  (OqIq  is  a  simple  frequency 
parameter  which  determines  the  analyzing  wavelet.  We  can  easily  see 
that  equation  (9)  satisfies  the  admissibility  condition  equation  (4).  The 
analyzing  bandwidth  of  the  wavelet  is  proportional  to  1/a,  thus  having  a 
constant  relative  bandwidth(BW),  that  is,  BW/(a)o/a)  =  constant.  This 

feature  is  also  reflected  in  the  narrow  time  window  at  higher  frequencyd.e., 
at  smaller  a  ).  In  general,  the  function  y/it)  is  selected  by  its  time  and 

frequency  localization  properties  [Daubechies,  1990]. 
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The  scale  change  may  be  performed  by  substituting  a  >=2'^  for  the 
computation  efficiency.  In  this  case  the  wavelet  family  is  V?  y/(2^  t), 
and  the  definition  of  the  wavelet  transform  becomes 

W(2i,6)  =  Jv(^)s(l)dt  (11) 

In  this  report,  the  unique  function  ^/(t)  can  be  selected  to  be  a  well-behaved 
modulated  Gaussian  function,  given  by  equation  (9). 
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III.  TIME-FREQUENCY  ANALYSIS 


Suppose  that  iff  is  any  basic  wavelet  such  that  both  Iff  and  its  Fourier 
transform  iff  are  window  functions  which  have  centers  and  half  widths  of 
time  and  frequency  domain  given  by  /*,  fi)*.  A/,  Aco  ,  respectively.  Then 
the  wavelet  transform  of  an  analog  signal  s(t)  is  given  as  follows. 

W(a,b)  =  J  s{t)  iffi^-^)  d/  (12) 

Equation  (12)  localizes  the  signal  with  a  time  window 
[b  +  at*  -  aAt,  h  +  a/*  +  aAi], 

where  the  center  of  the  window  is  at  b  + at*  and  the  width  is  given  by  2aAt. 
This  is  called  time  localization  in  signal  analysis. 

On  the  other  hand,  let  set  the  Fourier  transform  of  iff 

77(a))  =  IffiO)  +  0)*),  (13) 

where  Tf  is  the  Fourier  transform  of  iff,  then  Tj  is  also  a  window  function 
with  center  at  O*  and  width  by  Ao),  and  by  Parseval  identity,  the  wavelet 
transform  in  equation  (12)  becomes  [Chm,  1992] 

W(a,6)  =  ® -  {s((0)e’^^  Tf(aico-—)do)  (14) 

L  a 

•<•00 

where  the  phase  shift  of  is  determined  by  translation  along  time  axis. 
Equation  (14)  is  also  localized  information  of  spectrum  s(0))  of  the  signal 
s(t)  with  frequency  window 


Q) 

a 


1  .  fi) 
-Afi),  — 
a  a 


a 


where  the  center  of  window  is  at  0)  la  and  width  is  given  by  2£i(Oj  a. 
This  is  called  frequency  localization. 


From  equation  (12)  and  (14),  time-frequency  window  of  wavelet 
transform  is  as  follows. 

•  1  •  j 

[b  -I-  at*  -  flAr,  b  +  at*  aA/]x[— - Aft),  —  +  -Afij]  (15) 

a  a  a  a 


The  covered  area  by  time- frequency  window  is  the  multiplication  of  time 
and  frequency  bandwidth  arotmd  the  center  of  window.  The  time-frequency 
window  is  shown  figure  4.  Eventually  it  is  considered  positive  frequencies, 
i.e.,  a  >  0,  the  basic  wavelet  may  be  chosen  that  the  center  (0  of  is  a 
positive  number.  The  ratio  of  the  center  frequency  to  the  width  of  frequency 
band  is  given  by 


(0*  I  a  _  (0* 

21^0)  I  a  2A(0 


(16) 


which  is  independent  of  the  location  of  the  center  frequency.  This  is  called 
constant  percentage  band  width  analysis  or  constant-Q  frequency  analysis. 
In  this  report,  the  octave  band  is  used  for  the  analysis  of  vibration  signal. 
The  frequency  window  along  the  frequency  axis  narrows  for  large  center 
frequency  (O  !  a  and  widens  for  small  one  and  the  time  window  is  opposite 
to  frequency  window.  The  area  of  the  window  is  a  constant,  given  by 
4ArA(it}. 


12 


t 


bi  +  a/ 


l>2  +  ^2^ 


Figure  4.  Time-frequency  windows,  J]  >  02- 
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IV.  DISCKETE  WAVELET  TRANSFORM 


In  the  continuous  wavelet,  the  family  is  considered 

(17) 
a 

where  h  €  R,  a  €  Rf  with  a  ^  0,  and  is  admissible.  R  denotes  the  real 
line  and  R.f  denotes  the  positive  real  line.  It  is  considered  with  discrete 
values  for  a  and  b.  For  the  discretization  of  the  dilation  parameter,  we 
select  a  s:  Qq,  where  m  e  Z,  and  the  dilation  step  1  is  fixed.  For 
convenience,  it  assume  Oq  >  1.  If  Qq  =  1,  the  wavelet  transform  may  be 

similar  to  windowed  Fourier  transform.  For  m  =  0,  it  seems  natural  as 
well  to  dicretize  b  by  taking  only  the  integer  multiples  of  one  fixed 
where  b^  is  appropriately  chosen  so  that  the  V^(t  •  nb^)  cover  the  whole 
line.  We  arbitrarily  fix  6^)  >  0  in  this  report.  For  the  different  values  of  m, 
the  width  of  ^  times  the  width  of  y^(t)  defined  at  section 

III,  so  that  the  choice  b  -  nb^^aQ  will  ensure  that  the  discretized  wavelets 
at  level  m  cover  the  line  in  the  same  way  that  y^(t  •  nb())  do.  Thus  it  is 
selected  ,a  =  Oq,  b  -  where  m,n  range  over  Z,  and  Oq  >  1,  bQ  >  0 

are  fixed;  the  appropriate  choices  for  Oq  and  ^  depend  on  the  wavelet  yf 
and  the  characteristics  of  signal.  The  discrete  form  of  equation  (17)  is 
given  as  following  equation  [Daubechies,  1992] 

V'm.n(0  =  ¥("  ) 

% 

(18) 

=  yf{a^X  -  nbQ) 

For  the  computation  efficiency,  we  assume  that  0  =  2"’,  that  is,  Og  =  2, 
where  m  is  termed  the  octave  of  the  transform.  This  means  that  the 
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frequency  resolution  of  wavelet  has  a  octave  band.  The  integral  equation  (3) 
yields  a  wavelet  seri  3S  as  following  equation  by  using  equation  (18). 


At  the  discrete  wavelet  transform,  the  finite  energy  for  the  wavelet 
transform  is  not  equivalent  to  finite  energy  for  the  wavelet  series.  It 
depends  on  the  sampling  grid  as  well  as  the  function  In  addition,  it 

often  take  ^  to  be  a  multiple  of  a . 


.4^  - - 

W(2"',  n2"’)  =  -  n)s(t)dl  (20) 

A  logical  step  in  applying  the  theory  to  discrete  signal  is  to  discretize  the 
integral  in  equation  (20)  as  follows. 

W(2'",  02"')  =  ^  -  n)s(k)  (21) 

Octave  m  is  only  output  every  2"'  samples.  In  this  form  the  resulting 
algorithm  will  not  be  translation  invariant  [Mallat,  1982].  The  discrete 
wavelet  transform  is  highly  not  invariant  under  translations.  In  practice 
one  does  not  use  an  infinite  number  of  scales,  but  cuts  off  very  low  and  very 
high  frequencies. 
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V.  EXAMPLES  AND  DISCUSSIONS 


A  signature  generated  by  machinery  involves  many  informations  about 
its  operating  condition.  It  can  be  obtained  the  information  about  the 
operating  condition  of  machinery  by  applying  the  analysis  tools  for  the 
vibration  records.  Wavelet  transform  is  a  new  tool  that  is  particularly 
suited  for  time-frequency  analysis  of  nonstationary  or  stationary  signals. 
There  are  many  advantages  of  using  wavelet  transform  for  both  steady  and 
transient  signals.  We  will  discuss  the  performance  of  the  wavelet 
transform  by  using  simple  example  and  compare  with  pseudo  Wigner-Ville 
distribution(PWVD)  in  time-frequency  domain  analysis. 

Before  showing  some  examples,  it  is  necessary  to  discuss  how  best  to 
describe  the  results.  We  have  found  that  3-D  map  and  2-D  map  of  wavelet 
transform  are  a  useful  presentation  for  many  applications.  The  square  of 
the  amplitude  in  equation  (21)  has  the  unit  power/Hz.  The  distribution  of 
the  amplitude  square  over  the  individual  wavelet  and  position  can  now  be 
seen  as  figure  5.  If  the  length  of  sampled  data  is  shorter  than  the  wavelet 
size,  the  distribution  at  lower  wavelet  level,  i.e.,  low  frequency,  has  the 
value  at  only  one  position. 

In  order  to  describe  the  results  of  wavelet  transform,  we  use  the  3-D  and 
2-D(for  the  complicated  signal)  graphics.  And  for  graphic,  the  results  of 
that  is  distributed  and  reduced  to  256(3-D)  or  512(2-D)  data  point  along  the 
time  axis.  Frequency  axis  is  a  log  scale  (octave  scale  or  wavelet  level)  and 
time  axis  is  a  linear  scale  in  figures  of  wavelet  transform. 


A.  Harmozuc  wave  with  stepwise  fiiequeiicy  changes 

Figure  6  shows  (a)  the  pure  sine  wave  with  stepwise  frequency  changes 
500  Hz,  250  Hz  and  100  Hz,  (b)  its  PWVD  and  (c)  the  wavelet  transform.  The 
wavelet  transform  and  PWVD  well  represent  the  time  delay  and  the 
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frequency  components  of  signal.  The  wavelet  transform  is  a  result  with 
frequency  partitioned  into  consecutive  octaves  for  computational  effrciency. 
Then  the  magnitude  of  100  Hz  is  dispersed.  But  the  wavelet  transform 
clearly  describes  the  time  delay.  From  these  figures,  we  can  see  that  the 
PWVD  has  the  higher  resolution  than  the  wavelet  transform  in  frequency 
line.  However  the  computation  time  of  PWVD  is  extremely  higher  than  the 
wavelet  transform.  Figure  7  shows  the  wavelet  transform  of  the  sine  wave 
with  500  Hz  in  time  from  0.085  sec  to  0.17  sec.  the  wavelet  transform  well 
represents  the  time  delay  and  the  frequency  band  of  the  signal. 


Time 


Figure  5.  The  lattice  of  time-frequency  localization  of  wavelet  transform 

(N  =  Number  of  data  points) 


Magnitude 


&  Swept  harmonic  wave 


The  effect  of  sweep  rate  on  PWVD  was  well  investigated  at  Jeon  and 
Shin  (1993),  and  Shin,  Jeon  and  Spooner  (1993).  Figure  8.  shows  (a)  the 
swept  cosine  wave  with  the  sweep  rate  32  Hz/sec,  (b)  its  PWVD  and  (c)  its 
wavelet  transform.  From  this  figures,  we  can  see  that  PWVD  is  a  useful 
tool  for  analysis  of  the  signal  with  fast  frequency  change  in  time.  However 
if  the  records  of  signal  is  longer,  the  computation  time  may  be  much 
needed.  The  resiilt  of  the  wavelet  transform  does  not  clearly  represents  the 
sweep  condition  but  well  describes  the  change  of  frequency  range  with  time. 
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(a) 


Fig.8  (continued) 


(c) 

Figure  8.  Time-frequency  localization  of  the  PWVD  and  wavelet  transform: 
(a)  signal  8(t)  *  cos  {2%  32 1^),  (b)  its  PWVD  and  (c)  its  wavelet  transform  (fs  = 
256  Hz,  N  =  256). 
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C.  Harmonic  wave  with  itches 


The  interesting  phenomena  on  the  signal  with  an  abnormal  component 
as  a  fault  were  investigated.  Figure  9  shows  (a)  the  harmonic  wave  with 
glitch  at  a  small  region,  (b)  is  PWVD  and  (c)  its  wavelet  transform.  It  can 
be  seen  that  both  PWVD  and  wavelet  transform  of  the  signal  figure  9(a)  well 
represents  the  location  of  each  glitch  and  its  frequency  components. 
Figure  10  is  the  result  of  wavelet  transform  in  case  of  sampling  frequency 
8192  Hz.  From  these  figures,  it  can  be  clearly  seen  that  the  wavelet 
transform  is  very  useful  tool  in  the  analysis  of  the  signal  required  higher 
time  resolution.  Figure  9(c)  shows  the  glitch  components  more  clearly  than 
PWVD  because  the  wavelet  transform  has  the  very  narrow  time  window  in 
high  frequency  region.  This  characteristic  of  the  wavelet  transform  is 
useful  to  detect  the  fault  or  glitch  although  the  fault  is  small  and  to  monitor 
the  condition  on  any  vibrational  machinery  under  the  steady  operation 
condition.  Also  the  wavelet  transform  is  more  effective  for  the  analysis  of 
signal  which  the  time  record  length  is  long,  since  the  sweep  along  the 
frequency  line  is  octave  step  . 


Fig.9  (continued) 
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Figure  9.  Time-frequency  localization  of  PWVD  and  wavelet  transform:  (a) 
the  signal  sCt)  [Jeon  and  Shin,  1993],  (b)  its  PWVD  and  (c)  its  wavelet 
transform  (fg  »  256  Hz,  N  *  256). 


D.  Harmonic  wave  with  pulse 


Figure  11  illustrates  about  the  signal  including  small  changes 
generated  by  early  damage.  In  practice,  this  signal  is  not  given  by  this 
continuous  expression,  but  by  samples,  and  adding  a  5-function  is  then 
approximated  by  adding  a  constant  to  one  sample  only.  The  example 
signal  is  generated  by  following  equation. 

8(t)  =  sin(2x  100 1)  +  sin(2x  500 1)  +  1.5  5(t-0.049)  +  1.5  5(t-0.061)  (22) 

From  this  figure,  we  can  see  that  the  wavelet  transform  has  the  great 
advantage  for  the  analysis  of  the  signature  including  the  small 
disturbance.  PWVD  does  not  give  the  exact  information  about  time  delay  of 
pulse  but  well  represents  about  the  frequency  components  of  main 
signatures.  The  wavelet  transform  does  not  describe  the  magnitude  of  the 
main  frequency  components  but  well  represents  the  time  delay  of  pulse  and 
is  obtained  the  result  by  very  short  computation  time.  From  this  example, 
the  wavelet  transform  is  very  useful  tool  to  detect  the  fault  although  that  is 
very  small  region  on  time  axis. 


Fig.  11  (continued) 


28 


nude 


*<P 


(c) 

Figure  11.  Time-frequency  localization  of  PWVD  and  wavelet  transform: 
(a)  the  signal  s(t),  (b)  its  PWVD  and  (c)  its  wavelet  transform  (fs  =  8192  Hz, 
N=1024) 
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E.  Actual  pump  signal 


The  signal  of  actual  pump  was  measured  at  the  steady  state  speed 
without  and  with  abnormal  condition.  Figure  12  show  the  time  signal 
pattern  for  each  condition,  and  figures  13  and  14  show  the  results  of  wavelet 
transform  by  using  3-D  or  2-D  graphics.  Also  figure  15  are  the  results  of 
pseudo  Wigner-Ville  distribution.  These  figures  show  a  very  interesting 
results. 

From  these  figures,  we  can  easily  find  the  abnormal  condition.  At  13th 
frequency  band  of  figure  13(a)  and  (b),  the  results  of  wavelet  transform 
obviously  show  the  different  condition  and  pattern  for  each  case.  And  also 
at  8th  frequency  band,  the  small  changes  of  patterns  can  be  seen.  Figure 
14(a)  and  (b)  well  represent  the  operation  condition  although  the  shape  is 
more  complicate  and  we  can  see  the  abnormal  state  of  pump. 

In  this  case,  if  one  wavelet  to  have  13th  frequency  band  is  used,  it  gives 
good  results  for  diagnosis  and  vibration  condition  monitoring  as  very  short 
computation  time.  Figure  15  are  PWVD.  Also  PWVD  well  represents  the 
abnormal  condition  but  the  computation  time  is  very  long  than  that  of 
wavelet  transform.  In  the  case  of  wavelet  transform,  it  is  possible  to 
independently  compute  for  each  wavelet  level  different  from  PWVD.  The 
difference  of  computation  time  for  the  wavelet  transform  and  PWVD  is 
about  100  times  for  the  number  of  sample  data  N=2048  in  the  calculation  of  a 
whole  plane.  At  VAX3520,  the  computation  time  is  18  seconds  for  wavelet 
transform  and  about  30  minutes  for  PWVD. 

Especially  the  wavelet  transform  may  be  very  useful  to  detect  the  small 
disturbance  over  long  record  length  and  to  analyze  the  signals  which  has  a 
long  time  duration  and  intermittent  abnormal  condition  such  as  non¬ 
rotating  machinery.  The  advantage  of  the  Wigner-Ville  distribution  is  that, 
unlike  the  wavelet  transform  or  the  windowed  Fourier  transform,  it  does 
not  introduce  a  reference  function(such  as  wavelet  or  window  function  in 
Eq.(3)  and  (D)  against  which  the  signal  has  to  be  integrated  and  the  short 
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time  signal.  The  disadvantage  is  that  the  signal  enters  in  Wigner-Ville 
distribution  in  a  quadratic  rather  than  linear  way,  which  causes  many 
interference  phenomena  as  shown  in  references  Jeon  and  Shin  (1993),  and 
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(a) 

Fig.l3  (continued) 
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(b) 

Figure  13.  Time-frequency  localization  2-D  map  of  wavelet  transform  for 
Fig.  12(a)  and  (b),  respectively  (fs  =  10  kHz,  N=2048). 
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VI.  CONCLUSIONS 


The  wavelet  transform  has  been  investigated  and  applied  to  analyzing 
sampled  signal  and  the  actual  pump  signal.  The  results  of  this  research 
will  be  valuable  asset  for  the  analysis  of  vibration  records  and  condition 
monitoring  of  machinery.  The  following  conclusions  can  be  drawn: 

(1)  The  wavelet  transform  is  ideally  suited  for  portraying  the  wide-band 
transient  or  nonstationary  vibration  records  in  time*frequency  domain. 

(2)  The  wavelet  transform  has  a  great  advantage  to  detect  the  small 
disturbance  of  the  signal. 

(3)  The  computation  time  of  wavelet  transform  is  very  short  in  comparison 
with  other  time-frequency  localization  techniques. 

(4)  The  modified  Gaussian  wavelet  was  well  behavior  and  very  effective  to 
analyze  the  vibration  records. 

(5)  The  wavelet  transform  characterizes  the  time-frequency  localization  of 
the  signal  well  and  may  be  useful  tool  for  the  machinery  condition 
monitoring. 

(6)  If  the  wavelet  level,  that  is,  frequency  band,  is  selected  at  wavelet 
transform,  its  results  will  be  a  useful  tool  for  the  effective  pattern 
recognition  of  machinery  diagnosis  and  monitoring. 
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APPENDIX.  PROGRAM  LIST 


A.  Prognun  list  of  Wavelet  Traxisform(FORTRAN  77) 

c 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c  This  program  is  a  wavelet  transform  by  using  modified 
c  Gaussian  wavelet  (octave  sweep  or  1/3  octave  sweep), 

c 

c  Variables 

c 

c  n  =  number  of  input  data(n  must  be  the  power  of  2.) 

c  fs  =  sampling  frequency 

c  dt  s  sampling  time(time  interval) 

c  bw  s  cutoff  frequency  in  highpass  filter 

c  as  dilation 

c  b  s  translation 

c  w  =  coefficient  of  wavelet  transform  I  w(a,b)  I  **2 
c 

c  Array 

c 

c  x(i)  s  input  data  set 
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program  wavelet 


WAVELET  TRANSFORM 
by 

Dept,  of  Mechanical  Engineering 
Naval  Postgraduate  School 


c  xl(i)  s  filtered  input  data  set 
c  bk(i)  =  filter  weight 
c 
c 

dimension  x(8192),bk(4100)ptl(8192) 
complex  ai.siun 
character*  16  inname, outname 
character  as 
c 

ai=cmplx(0.,l.} 

pi=atan(l.)*4. 

c 

print*, ’===================:=============! 

print* 

print*,'  WAVELET  TRANSFORM' 

print* 

print*, 's=:===!=:=sss=!s:===s:=ss:=a:====s=s=!===!— 
print* 

print*,'What  is  an  input  filename  ?’ 
read(*,10)  inname 

print*, "What  is  an  output  filename  ?' 
read(*,10)  outname 
10  format(al6) 

c 

c  select  the  sweep  method 

c 

print*,'What  do  you  want  the  sweep  method  ?' 
print*,'  if  1/1  octave,  input  1' 
print*,'  if  1/3  octave,  input  2' 
read(*,*)  mswp 

c 

c  read  the  input  data  file 

c 

call  indata(inname,x,fs,n,nn) 
c 
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dtsl./fs 


c 

call  smean(n,x} 
c 

print*,'Do  you  want  to  apply  highpass  digital' 
print*, 'filter  to  the  original  data  ?  (y/n)' 
read(*,20)  as 
20  format(al) 

if  (as.eq.'Y'.or.as.eq.'y’)  then 

print*, 'Enter  the  cutoff  frequency  of 
print*,'the  digital  highpass  filter  (in  Hz)’ 
read(*,*)  bw 
endif 
c 

c  signal  modifications 

c 

c  Application  of  highpass  digital  filter 

c 

if  (as.eq.'Y’.or.as.eq.'y')  then 
mo=n/2 

c  calculate  the  filter  weighting 

call  filter(mo,bw,dt,bk) 

write(*,*)  'finished  filtering' 

c  pass  the  highpass  filter 

do  160  isl^ 

100  xl(i)=0. 

(k)  200  isl,n 

do  170  ks.mo,mo 
i=k 
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if  (k.lt.O)  then 
j=k*(-l) 

endif 

j*j+l 

ll»i-k 

if  (ll.lt.  1)  then 
ll=ll-t-n 

elseif  (ll.gt.n)  then 
llsll-n 

endif 
bb=-bk(j) 
if  (k.eq.O)  then 
bb=l-bk(j) 

endif 

xl(i)=xl(i)+bb*x(ll) 

170  continue 

200  continue 

do  210  isl,n 
210  x(i)=xl(i) 

endif 

if  (mswp.eq.l)  then 
nk=nn 
facssl. 

elseif  (mswp.eq.2)  then 
nkss3*nn 
facs3. 
endif 

nrsnk-int(fac)-«- 1 

open(7,f!lesoutname,status='new') 

c 

c  calculate  wavelet  transform 
c 

tini=0. 
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ttimesdt*!! 
fDss2.''‘pi 
write(7,*)  tini 
write{7,*)  ttime 
write(7,*)  nn,index 

do  1000  isl,nr 

fio=2**((nk-i+l)/fac) 

a=fio*dt 

nm=int(ttime/a) 
if  (nm.eq.O)  nmsl 
ddtsttime/nm 
sr=sqrt(-alog(0.0 1  )*2)*fio 
max=int(sr) 

if  (max.ge.2’'‘n)  maxs:2*n 

do  900  j=l,nm 

suin=cmplx(0.,0.) 

bsddt/2.+ddt*float(j-l) 

jk=int(b/dt)+l 

do  800  ks-max,max 

t=k*dt/a 
klsjk+k 
if  (kl.lt.  1)  then 

if  (kl.lt.-n)  then 
kla;kl+2*n 

else 

klskl+n 

endif 

endif 

if  (kl.gt.n)  then 

if  (kl.ge.2*n)  then 
klsskl-2*n 
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else 


klskl-n 

endif 

endif 

8uin=s\im+x(kl)''‘cexp(-ai*f0’*‘t)*exp(-t*t/2.)*dt 

800  continue 

w=cabs(sum)/sqrt(a) 
w=abs(w)*abs(w) 
write(7,*)  ij.w 

900  continue 

1000  continue 
cIose(7) 
stop 
end 
c 
c 
c 

c  SUBROUTINES 
c 

C - 

subroutine  indata(inname,x,fs,n,nn) 

dimension  x(*) 
character*  16  inname 
c 

open(5,file=inname,status='o]d') 
read(5,*)  fs,n 
do  100  j=l4i 

read(5,*)  t,x(j) 

100  continue 
close(5) 
ft*fs/2. 
do  is  1,20 
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ii*2**i 
disfloatdi) 
if  (di.ge.ft)  then 
nnsi 
go  to  200 

endif 
end  do 

200  continue 

return 
end 
c 
c 

subroutine  smean(n,x) 
c 

c  This  subroutine  calculates  and  removes  the  mean  value  of 

c  the  sampled  data, 

c 

dimension  x(*) 
real  meanv 
c 

asumsO. 
do  100  i=l,n 

asum=asum+x(i) 

100  continue 

meanvsasum/n 
do  200  i=l,n 

x(i)=x(i)-meanv 
200  continue 

return 

end 

c 

c 

subroutine  filter(mo,b,t,bk) 
c 

c  Routine  generates  FIR  filter  weights. 

48 


M<jthod  devised  by  Potter,  Bickford  and  Glaze. 

There  are  a  total  of  2M+1  weights...filter  generates  M+1. 

—  variables  — 

t  s  the  sampling  interval  in  second, 
b  =  cutoflKhalf-power  point)  of  the  filter  in  Hz; 
must  be  on  the  range  firom  o  to  l/2t. 

Results  are  stored  in  bk 

•  Note;  in  the  case  of  highpass  filter,  the  value  of 
weight  bO  miist  use  l*bO  instead  of  bO. 

dimension  bk(*),d(3) 

datad0/0.35577019/,d(iy0.2436983/,d(2)/0.07211497/, 

*  d(3)/0.00630165/ 
pi=atan(l.)*4. 
msmo 

first  generate  plain  boxcar  weights 

facts2.*b*t 

bk(l)ssfact 

fact=fact*pi 

do  5  isl,m 

fi=i 

bk(i+l)=sin(fact*fi)/(pi*fi} 
trapezoidal  weighting  at  end 
bk(m+ 1  )=bk(m+ 1  )/2 . 

Now  apply  the  Potter  p310  window 

sumgsbkd) 

do  15  isl,m 

sum=dO 

fact=pi*float(i)/float(m) 
do  10  k=l,3 

8um=sum+2.*d(k)*cos(fact''‘float(k)) 

bk(i+l)=bk(i+l)’'‘sum 

sumgs=8umg+2.*bk(i+l) 
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ml=m+l 
do  20  isl,ml 
20  bk(i)=bk(i)/siimg 
return 
end 
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&  Pmogram  List  of  wavelet  transform{MA'nLAB) 


% 

% 

Wavelet  Transform 

% 

% 

by 

% 

Dept,  of  Mechanical  Engineering 

% 

Naval  Postgraduate  School 

% 

% 

Ver.  1.0  Aug.  5  1993 

% 

% 

By  using  MATLAB  Ver.  4.0 

% 

% 

“  Variables 

% 

infile 

=  input  filename 

% 

fs 

=  sampling  frequency 

% 

n 

=  the  number  of  sampled  dat 

% 

( n  must  be  n-th  power  of  2.) 

% 

avg 

s  mean  value  of  sampled  data 

% 

fact 

=  the  factor  of  sweep  rate 

% 

x(i) 

=  the  magnitude  of  signal 

% 

z(ij) 

=  the  result  of  wavelet  transform 

% 

load  c:\users\infile 

f8=infile(l,l); 

n=iniile(l,2); 

dt=l/fs; 

x=infi1e(2;n-fl,2); 

asumsO; 

% 

%  remove  mean  value 

% 

for  jl=l:n 

asum=asum+x(jl); 
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end 


avgsasuzn/n; 
for  ils:l;n 

x(il)=:x(il)-avg; 

end 

% 

%  Decide  the  sweep  rate 
% 

factsinputCsweep  method:  1/1  octave  *1,  1/3  octave  =  3  '); 

% 

ft=fs/2; 

forilssl:20 

if  ii  >s=  ft 
nn=il; 
break 
end 
end 

if  fact  ss  1 
nksnn; 

else 

nk=:fact*nn; 

end 

nrs:nk>fix(fact)+ 1 ; 

% 

%  calculate  wavelet  transform 

% 

tini=0; 

ttime=:n*dt; 

fDs2*pi; 

jmasfix(ttime/(2*dt)); 
wt=zeros{nr  jm); 

for  ila  1 :  nr 

ios2^((nk-il+l)/fact); 

asio^dt; 

nm=fix(ttime/a); 
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if  nmssO 
nm  s  1; 
end 

ddtsttime/nm; 
sr=8qrt(-log(0.01)'''2)*io; 
inxsfix(sr); 
if  mx  >s  (2*n) 
mx=2*n; 
end 

for  jsl:nm 
8um=0+i*0; 
b=ddt/2+ddt*(j-l); 
jk=fix(b/dt)+l; 
for  k=-mx:mx 
t=k*dt/a; 
kl=jk+k; 
if  kl<l 
if  kl  <  -n 
kl=kl+2*n; 
else 

kl=kl-fn; 

end 

end 

if  kl>n 

if  kl  >=  2*n 
kl=kl-2*n; 
else 

kl=kl-n; 

end 

end 

xl*0.-i*fDn; 

8uni*8um+x(kl)*exp(xl)*exp(-t*t/2)*dt; 

end 

w®abs(8um)/8qrt(a); 

ws:ab8(w)*abs(w); 
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wt(ilj)=ww; 


% 

%  Plotting 
% 

lines:  128; 

zszero8(nr,line); 

foril»l:nr 

ioss2'^((nk-i+iyfact); 

assio*dt; 

k=lix(ttiine/a); 

ifk==0 

ksl; 

end 

if  k  <-  line 
stepsline/k; 
k2=K); 
for  ii=l:k 
kl=k2+l; 
k2a:fix(8tep*ii); 
st=!k2*kl+l; 
if  8t<s:8tep 
cfsstep/st; 
else 

cf=st/step; 

end 

for  ijskl:k2 

z(i  1  ,ij  )s  wt(i  1  ,ii)/cf7  st; 
end 
end 
else 

step=k/line; 

k2=>0; 

for  ii=l:line 
kl=k2+l; 
k2=!fix(step*ii); 
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st*k2-kl+l; 

ssumsO; 

for^skl:k2 

ssum=ssum+wt(il,ij); 

end 

if  St  <a  step 

ssumsssuni*8tep/st; 

else 

ssumsssum^st/step; 

end 

z(il,ii)=:BSum; 

end 

end 

end 

xmasmax(z); 
xp=max(xma)*1.2; 
dtl=ttime/(line-l); 
for  ksliline 

xvalues:(k-l)*dtl; 

xx(k)=xvalue; 

end 

for  k=l:nr 
yy(k)=k; 
end 

surftxx,yy,z) 
xlabeK'Time  (sec)’) 
ylabelCFrequency  step  ') 
zlabel(’Amplitude') 
axis([tini  ttime  1  nr  0  xp]) 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


20 

c 

c 

c 


C.  Plotting  Pkx»gram  List  (3  Dimension) 

•  FOTRAN  77  and  CA-DISSPLA  Graphic  Package  • 

program  plot 

This  program  uses  the  graphic  package  CA-DISSPLA  to 
plot  the  results  of  wavelet  transform 


tmax 

=  time  record  length 

tint 

=  initial  time 

fmin 

s  start  frequency  step 

fmax 

s  stop  firequency  step 

nn 

=  the  number  of  the  frequency  step 

fac 

=  index  of  sweep  method 

fname 

=  input  data  filename  generated  by  main  source 
program 

dimension  rr(32768),wt(50,256) 
character’*'25  fname 

write(*,*)  'input  file  name  ?’ 
read(*,20)  fname 
format(a25) 

data  distribution  or  reduction  for  3-D  graphics 
lines256 


open(8,filesfname,statuss'old') 

read(8,’'')  tini 
read(8,*)  tmax 
read(8,*)  nn.n.fac 
nk=nn*int(fac) 
nrssnk-int(fac)+l 
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ddtstmax/n 
do  1000  i=l4ir 

fio=2.**((nk-i+l)/fac) 

a*fio*ddt 

k=int(tmax/a) 
if(k.eq.O)  kssl 
if  (k.le.line)  then 

step=:float(line)/float(k) 

k2s0 

do  950  ii=l,k 

read(8,*)  il  jl,w 
kl=k2+l 
k2=int(step*ii) 
st=float(k2-kl)+l. 
if  (st.le.step)  then 
coe=step/st 
else 

coesst/step 

endif 

do930ij-kl,k2 
wt(i,ij)=w/coe/st 
930  continue 

960  continue 

else 

step=float(k)/float(line) 

k2s0 

do  800  iislJUne 
kl=k2+l 
k2=int(step*ii) 
sumsO. 

8t=float(k2>kl>+l . 
do  750  \j=kljc2 
read(8,*)  il  jl,w 
sumssum+w 
750  continue 
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if  (st.le.step)  then 
8umssum'*‘step/st 
else 

sum«sum*st/step 

endif 

Mrt(i,ii)»sum 


800 

continue 

endif 

1000 

continue 

close(8) 

smaz=0. 

do  2000  i=l,nr 

do  2000  jsljine 

if  (smax.le.wt(ij))  smax=wt(ij) 
2000  continue 

write(*,*)  'maximunis'jsmax 
write(*,*)  'input  the  maximum  scale  ?' 
read(*,*)  fact 

do  120  isl,nr 
do  100  j=l, line 
ks:line*(i-l)+j 
rr(k)=wt(ij) 

100  continue 

120  conunue 

ddsmod(nr,2) 
if  (dd.eq.O)  then 
nrssnr+l 

kk=(nr-l)*line+l 
klsnr^line 
do  150  ip=:kk,kl 
150  rr(ip)=!0. 

endif 


58 


c  plotting 
c 

call  pdev('ln03',ieer) 
call  hwshd 
call  swissm 

call  shdchrOO., 1,0.002,1) 
call  heighUO.lS) 
call  physor(l.l,1.2) 
call  area2d(5.5,6.75) 

call  messagCWAVELET  TRANSFORM  $’,100,1.1,8.2) 
call  blsiu* 

call  volm3d(8.,8.,9.) 

call  xSnameCTime  (sec)  $',  100) 

call  y3naine('Frequency  step’,  100) 

call  z3name(’Amplitude  $',  100) 

call  zaxangOO.) 

finaxafloat(nr) 

f8tepss(fmax-l)/4. 

call  graf3d(tini,tmax/4.,tmax,l.,fstep,fmax,0.,’SCALE’,fact) 

call  suxmat(rr,l,line,l,nr,l) 

call  end3gr(0) 

call  endpKO) 

call  donepl 

end 
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D.  Plotting  Program  List  (2  Dimension) 

•  FOTRAN  77  and  CA-DISSPLA  Graphic  Package  • 

c 

program  levelplot 
c 

c  This  program  uses  the  graphic  package  CA-DISSPLA  to 
c  plot  the  results  of  wavelet  transform  with  respect  to  each  level, 
c 

c  tmax  s  time  record  length 

c  tint  =  initial  time 

c  fmin  =  start  frequency  step 

c  fmax  =  stop  frequency  step 

c  nn  =  the  number  of  the  frequency  step 

c  fac  s  index  of  sweep  method 

c  inname  =  input  data  filename  generated  by  main  source 
c  program 

c 

dimension  x(512),y(512),wt( 50,256) 
character’"  25  inname 
c 

writeC",’")  ‘input  file  name  ?' 
read(’",20)  inname 
20  format(a25) 

c 

c  data  distribution  or  reduction  for  2-D  graphics 

c 

line=512 

open(8,file=fname,status='old') 

read(8,'")  tini 
read(8,’")  tmax 
read(8,'")  nn,n,fac 
nks:nn'"int(fac) 
nr=nk-int(fac)+l 
ddt=tmax/n 
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1000  issl^ 
fios2**((nk-i+l)/fac) 
a=fio*ddt 
k=int(tmax/a) 
if  (k,eq.O)  k*! 
if  (k.le.line)  then 

8tepsfloat(lineVfloat(k) 

k2=0 

do  950  ii=l 
read(8,*)  iljl.w 
kl=k2+l 
k2=int(8tep*ii) 
st=float(k2-kl)+l. 
if  (st.le.step)  then 
coesstep/st 
else 

coesst/step 

endif 

do  930  ijBkl,k2 
wt(i,ij)*w/coe/st 
continue 
continue 
else 

step=float(k)/float(line) 

k2^ 

do  800  ii=l4ine 
kl=k2+l 
k2=int(step'*ii) 
sum=0. 

8t=float(k2-kl)+l. 
do  750  y=kl  4(2 
read(8,*)  iljl,w 
8uni=suin+w 
continue 
if  (stJe.step)  then 
8um=sum*step/8t 
6  1 


else 

sumssum  "‘st/step 
endif 

wt(i,ii)=sum 


800 

continue 

endif 

1000 

continue 

close(8) 

smaxsO. 

do  2000  i=l,nr 

do  2000  jsl, line 

if  (smax.le.wt(ij))  smaxswtCiJ) 
2000  continue 

write(*,*)  'maximum=',smax 
write(*,*)  'input  the  maximum  scale 
read(*,*)  smax 
c 

c  normalizing 

c 

do  2500  i=l,nr 
do  2500  jad.line 

^00  wt(ij)=wt(ij)/smax 

c 

c  plotting 
c 

call  pdev(’ln03',ieer) 
call  hwshd 
call  swissm 

call  shdchrOO., 1,0.002,1) 
call  height(0.15) 
call  page(8.5,ll.) 
call  physor(1.5,1.5) 
call  area2d(4.8,4.7) 
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2700 


3000 


call  xnameCTime  (sec)  $',  100) 

call  yname( 'Frequency  step’,  100) 

call  headinCWAVELET  TRANSFORM  $’,100,1.1,1) 

call  thkfrm(O.Ol) 

call  yaxangCO.) 

call  graf[tini,tmax/4.,tmax,l.,l,,nr+l) 
call  grid(l,l) 
do  3000  issl,nr 
do  2700  j=l,line 
x(j)s:dt*float(j) 
y(j  )=float(i)+wt(i  j) 
continue 

call  curve(x,y,line,0) 

continue 

call  endpKO) 
call  donepl 
stop 
end 
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